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We investigate the possibility to assist the numerically ill-posed calculation of spectral properties 
of interacting quantum systems in thermal equilibrium by extending the imaginary-time simulation 
to a finite Schwinger-Keldysh contour. The effect of this extension is tested within the standard 
Maximum Entropy approach to analytic continuation. We find that the inclusion of real-time data 
improves the resolution of structures at high energy, while the imaginary-time data are needed 
to correctly reproduce low-frequency features such as quasi-particle peaks. As a nonequilibrium 
application, we consider the calculation of time-dependent spectral functions from retarded Green 
function data on a finite time-interval, and compare the Maximum Entropy approach to direct 
Fourier transformation and a method based on Pade approximants. 
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I. INTRODUCTION 

A common strategy to circumvent the oscillatory con- 
vergence of integrals in interacting quantum field the- 
ories is the Wick rotation to imaginary time. Apart 
from possible fermionic sign problems, quantum Monte 
Carlo (QMC) simulations are well-conditioned on the 
imaginary time axis, and simulations at non-zero tem- 
perature using the Matsubara technique are straight- 
forward. However, extracting dynamic properties from 
imaginary-time data of the single-particle Green's func- 
tion G{— it) = —(T t c{t)c^{Q)) is an ill-posed problem: 
the singular integral equation 

iG(-ir) = / d^H^^ (1) 

has to be solved for the spectral function A(ui). 

If the real-time data for the retarded Green's function 

G ret (i,t') = -i6(t-t')({rf(f),rf t (*')}) (2) 

were known on the entire real-time axis, the spectral 
function could be obtained from a simple inverse Fourier 
transform (in equilibrium, G lct (t, t') depends only on the 
time difference t — t'). The Fourier transform is not 
ill-posed, due to unitarity. Unfortunately, calculating 
these real-time data for an interacting system is difficult. 
Monte Carlo techniques cannot reach long times, due to 
a dynamical sign problem which grows exponentially as 
a function of the maximum real time i max - Nevertheless, 
as shown in Refs. [IHS] the real-time Green function up to 
some finite time t max can be computed accurately using 
continuous-time Monte Carlo algorithms [1HZ] , or, in cer- 
tain parameter regimes, using perturbative weak- [8] or 
strong-coupling methods [9j [10] . This raises the question 
if the information contained in the real-time correlators 



for t, t' < i max can be exploited to obtain reliable spectra 
using a suitably adapted analytical continuation proce- 
dure. 

A real need for the "analytical continuation" of real- 
time Green functions arises in the study of nonequilib- 
rium properties of correlated systems, e. g. lattice sys- 
tems perturbed by a quench or some field pulse [TlTfPB"] . 
In order to characterize the relaxation of these systems, 
it is often useful to define a time-dependent spectral func- 
tion 
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A(u,t) = Im / dt'e lu(t -^G Ict {t',t), (3) 

*" Jt 

where G ret (t',t) now depends on both times individually 
due to the loss of time-translation invariance. A(u>,t) 
is not assured to be positive, but it typically becomes 
positive a short time after the perturbation. (In par- 
ticular, A(uj,t) is positive for any quasi-stationary state. 
When A(uj, t) is constant over a time window of width 
At, then A(ui,t) must be positive after averaging over 
a frequency window of Alo cx 1 /At.) Under certain as- 
sumptions, A(u>,t) is related to the photoemission and 
inverse photoemission signal [14L 115] . Furthermore, in an 
equilibrium system, Eq. ^ reduces to the familiar spec- 
tral function defined in Eq. (JlJ. The challenges in the 
evaluation of A(ui,t) are the same as those for equilib- 
rium spectra mentioned above: In practice, simulation 
results will be limited in time, so direct Fourier trans- 
formation will lead to oscillations or a smearing-out of 
spectral features. Therefore, a second question which we 
want to address is whether Maximum Entropy or Pade 
approaches can be used to improve the quality of these 
time-dependent spectra. 

The paper is structured as follows: in Sec.[n]we present 
and test the Maximum Entropy approach for Keldysh 
Green's functions. In particular, Sec. |IIB| investigates 
the singular values of the kernel matrix and their depen- 
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FIG. 1: (color online) The Schwinger-Keldysh contour C. Ar- 
rows denote the direction of time-ordering <c- 



dence on the choice of data points, in Sec. |II C| we test 
the analytical continuation procedure for exactly known 
Green's functions, and in Sec. |IID| we compare the Maxi- 
mum Entropy result to spectra obtained by Fourier trans- 
formation. A generalization of the Pade analytical con- 
tinuation procedure to nonequilibrium Green functions is 
introduced in Sec. |III| In Sec. |IV| we apply the Maximum 
Entropy method to real-time simulation data of equi- 
librium and nonequilibrium systems, and compare the 
method to direct Fourier transformation and the Pade 
procedure. Sec. [V] gives a summary and conclusion. 



II. MAXIMUM ENTROPY APPROACH 

A. Formulation of the problem 

In practice, real-time simulations are carried out on 
the Keldysh contour C [To] IT?] illustrated in Fig. [lj The 
Green's function is a function of two time indices t and 
t', which may lie on the upper real, the lower real, or the 
Matsubara branch of the contour. Using the Keldysh- 
contour time-ordering <c, the contour-ordered equilib- 
rium Green's function is related to the spectral function 
through 



iG(z,z') = (T c d(z)d\z')) 



with the contour-ordered Fermi factor 




if z' < c z, 
else, 



(4) 



(5) 



and the times z, z' located on C. Obviously, Eq. Q is a 
generalization of Eq. |TJ), and extracting A(w) from any 
data set G(zi,Zj) is still an ill-posed problem when the 
times Zi are restricted to the Matsubara branch, or to 
times smaller than some finite t max - 

The Maximum Entropy analytic continuation proce- 
dure involves the inference of a most probable spectral 
function A(uj) among those which are compatible with 
the simulated data. The problem can be expressed in 
the matrix form 



where D represents simulation data, if is a linear opera- 
tor, given by the matrix kernel, applied to the unknown 
spectral function A. The finite set D can be interpreted 
as a snapshot from a Gaussian random variable with co- 
variance C, in the case of a Monte Carlo simulation |19j . 

The essential idea for inferring a most probable A from 
D has been outlined in Ref. [15]. An entropy functional 
aS'LA] is added in order to regularize the minimization of 
X 2 [-D,^4] with respect to A, where 
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X 2 [D, A] :— - (KA - D) T C-\KA - D) 



(7) 



The most informative scalar a can be determined by 
a systematic application of Bayesian inference (Classic 
MaxEnt), or can be averaged over the corresponding 
Bayesian probability distribution (Bryan's MaxEnt). 

The standard numerical algorithm for the maximum 
entropy inference of spectral functions was developed by 
Bryan |21] . The kernel is decomposed via a singular- value 
decomposition 



K = VZU 1 



(8) 



where £ = diag (ci, 02, . . . ), o\ > 02 > ■ • • > 0. 
Through the matrix products, each singular value Oi is 
associated with a direction in ^4-space and a direction 
in D-space. The former is given by one of the "basis 
functions" for the spectral function and the latter is rep- 
resented by an entry of V. 

If <jj is large, it provides a channel which transports 
a comparably large amount of information about A. If 
crj is small, not much information can be gained for the 
corresponding direction in A-space, and the Bayesian ap- 
proach will not modify a default spectrum with respect 
to that direction, due to a lack of evidence. A small value 
of a can only be compensated by small error bars for the 
corresponding D-direction. Hence, the shape of the sin- 
gular value distribution is an important indicator for the 
structure of the inverse problem. 



B. Singular value distributions 

In this section we analyze the singular value structure 
of the integral kernel K in Eq. Q , for various data sets D 
(the set A is always given by some appropriate discretiza- 
tion of the spectral function). In panel (a) of Fig. [2] the 
singular value distribution is plotted for the usual infer- 
ence from imaginary-time data (solid lines), for inverse 
temperature (3 = 10. In this case, the data set D = Amag 
consists of the values G(—iTj,0) on N equidistant points 
on the imaginary time branch, 



D 



N,/3 



{G(—iTj,0) I Tj = 0j/N,j = 0, . . . , AT — 1}. (9) 



D = K ■ A, 



(6) 



The singular values are seen to decay exponentially. 

Next we define a suitable data set for extracting A(ui) 
from pure real-time data. We note that the imaginary- 
time set of input data contains real-time information for 
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FIG. 2: (color online) Singular value distributions for f3 = 10 
and indicated number of time-points N. Panel (a): Mat- 
subara and real-time Green's functions. Panel (b): Keldysh- 
offdiagonal Green's functions. 



both positive (forward direction) and negative (back- 
ward direction) frequencies. In the case of real-time 
data, the situation is different: If only one time order- 
ing of d, is taken into account, the Maximum Entropy 
method only yields information about the positive or neg- 
ative frequency range, i. e. the occupied and unoccupied 
part of the spectrum. It is thus important to consider 
"lesser" and "greater" (or retarded) Green's functions. 
We choose (z, z') values that correspond to both elec- 
trons propagating forward in time starting from if = 0, 
G > (t, 0) = — (0)) , and backward in time start- 

ing from t' = i max along the upper real-time branch, 
G < (i,t max ) = i(d* (t max )d(t)) , using an equidistant grid 
of annihilation times, 

= {ReG> (tj , 0) , ImG> {t 3 , 0) , ReG< (tj , t max ) , 

ImG< {tj, t max ) | tj = jt max 4/N, j = 0,...,N/4- 1}. 

(10) 




singular value index 



FIG. 3: (color online) Dependence of real-time singular value 
distributions on the number of input data iV in £ , ro ^ 1 max . 

Due to translational invariance, the d^ operator need not 
be shifted in time. Since we use both the real and imag- 
inary parts of G > (t, 0) and G < (t, £ max ) each time step 
contributes two real variables to the inference process, 
and we always use N for the the total number of real 
variables below. 

In contrast to the case of imaginary-time input, the 
kernel for pure real-time propagators has a distinguished 
plateau in its singular value distribution, even for rather 
small < max = 2 (dashed line in Fig. [2]). As £ max is in- 
creased (i max = 10, dotted line), the plateau broadens, 
and more equally dominant A-directions appear. This is 
because in the limit i max — > oo the (no longer ill-posed) 
unitary limit <7j = 1 is approached. 

We next consider singular-value distributions for 
Keldysh-offdiagonal propagators. For symmetry reasons 
we can restrict ourselves to mixed Green functions from 
the upper Keldysh contour to the imaginary branch, 
G^(r, t cnt ). We consider hxed values of i cut . t cut = is 
equivalent to the usual kernel for the Matsubara Green's 
function. As t cut is increased, the singular value dis- 
tribution is broadened (panel (b) of Fig. [2| , but not as 
dramatically as in panel (a) of Fig. [2] Keeping t = r cut 
fixed yields a structure more similar to the real-time dis- 
tributions in panel (a) (r cut = is exactly the same). 
However, it appears that the plateaus in panel (a) can- 
not be exceeded. 

Figure [3] shows the dependence of the real-time 
singular-value distribution on the number N of time- 
points, for a Keldysh contour length i max = 2. We 
find that at low singular value indices, apart from a fi- 
nite offset, the same shape is obtained. In particular, 
the same number of singular values is associated with 
the plateau. Furthermore, the initial descent from the 
plateau is also identical. However, the singular value dis- 
tribution for N = 100 continues to decrease smoothly, 
whereas the corresponding curve for N — 20 abruptly 
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jumps to the lowest singular value which is of the order 
of 10~ 15 . A similar behaviour is also found for other 
correlators within the Keldysh contour. 

We conclude that increasing N does not extend the 
width of the plateau (once the latter has been estab- 
lished), but merely adds further singular values to the 
rapidly decreasing tail. Since this decrease is exponen- 
tial, an increase in N is only useful if the added singular 
values are larger than the order of magnitude e of the 
numerical error of the data. This dependence is indeed 
observed in the data analyis presented in the following 
subsections. 



C. Tests of equilibrium spectral functions 

We will first analyze the behaviour of the continuation 
procedure using artificial data sets corresponding to given 
spectral functions with sharp peaks or band edges. For 
this purpose, uncorrelated data sets taken from the exact 
contour Green function are studied, and we assume a 
uniform error distribution e, i.e., the covariance matrix is 
taken to be of the form C — diag (e 2 ). The performance 
of the Maximum Entropy procedures on more realistic 
data sets will be investigated in Sec. |TV} 



0.3 



1. Rectangular spectrum 

As a first test case we consider the rectangular spectral 
function 



1, if|w|<2, 
0, else. 



(11) 



It can be expected that the sharp band edges will be dif- 
ficult to infer from any finite data set. It is well-known 
that even the inverse Fourier transform, i.e. the unitary 
limit i max — > oo converges only slowly, and that the con- 
vergence is not point-wise. Hence, we consider this an 
interesting test case for the Maximum Entropy method. 
The analytic continuation procedure is tested for inverse 
temperature /3 — 10 and for real-branch lengths i max = 2 
and t max = 10. The fake variance of all C-contour corre- 
lator estimators is set to e 2 = 6 • 10~ 14 . 

Figure [4] compares the exact spectral function to the 
A{uS) obtained from the analytical continuation for dif- 
ferent data sets ^ and (10). The total number of real 
variables is kept constant at N = 100, i.e., we use either 
100 equidistant time steps on the imaginary axis, or 25 
time steps in the real-time set -D^ai" 1 "*- (We found that 
the restriction of real input data to the real or imaginary 
parts of G(t, t') does not contain sufficient information 
for an analytical continuation.) 

In Fig. [4] it can be seen that using a broad Gaussian 
default model, the width of which has no significant in- 
fluence on the results, the maximum entropy solution 
shows Gibbs ringing artifacts. The frequency of these 
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FIG. 4: (color online) Comparison of spectra at /3 = 10 for 
N — 100 data points located on the Matsubara contour, or 
on the upper Keldysh contour up to time i ma x = 2, 10. 



oscillations is a measure of the accuracy of the inferred 
spectrum. Surprisingly, already the real-time data from 
a short Keldysh branch, t max — 2, yield more accurate 
results than the j3 = 10 Matsubara branch data. As £ max 
is increased, the approximation of the spectrum becomes 
even better. 

Further increasing the number of data points N does 
not significantly change the above results. However, low- 
ering the variance e 2 and then raising N systematically 
yields more accurate spectra for both, Matsubara and 
Keldysh data, because in this case the requirements for 
the identity theorem of complex analysis which guaran- 
tees uniqueness arc approached systematically. However, 
the convergence is always exponentially slow, as pointed 
out in the discussion of the dependence of the singular- 
value distributions on TV. 



2. Asymmetric triangle 



We next investigate 



AriagM 




I- 



if \UJ 

else 



21 < 1, 



(12) 



as an example for spectra which are not particle-hole 
symmetric and for which the convergence of the inverse 
Fourier transform is more rapid than for the rectangu- 
lar case. The variance is again set to e 2 = 6 • 10 -14 . 
Results for this scenario are shown in Fig. [5j Judged by 
the peak position and the steepness at the discontinuity, 
even the t max = 0.2 Keldysh spectrum is already slightly 
better than the Matsubara spectrum. This indicates that 
real-time data are good for resolving high-energy features, 
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FIG. 5: (color online) Comparison for the triangular spec- 
trum. 

since the triangle is shifted to relatively high energies. 
As in the case of the rectangular shape, the spectrum 
improves as t max is increased. 



3. Multiple peaks 

As a next step towards more realistic situations we turn 
to a spectrum with a sharp resonance at uj = and two 
side bands. We model this by superimposing Gaussians, 

-4pcakM := X] c sb&r sb ,anM + C r cs5<r r c S ,o(w), (13) 
a = ±l 

with c ros = 0.1, c sb = 0.45, £1 = 2.0, a sh = 0.5, a Tes = 

0.05; 9a,x(x) - 7=f exp 

As one can see in Fig. |6j real-time data on relatively 
short contours (£ max = 2, 10) are not particularly use- 
ful, when sharp low-frequency features such as the given 
"quasi-particle peak" should be extracted. The Mat- 
subara data and the short-£ max real-time data yield side 
bands of similar quality, while the imaginary-time data 
produce a comparable or even better reconstruction of 
the sharp resonance. Increasing i max mainly improves 
the resolution of smooth high-energy features (t max = 10, 
N = 400). We also observe that using n — N/A = 5 
equidistant time steps on the £ max = 2 real contour 
yields a better spectrum than a similar data set on the 
t max = 10 contour. This finding is compatible with the 
earlier suggestion that as a function of N, a full establish- 
ment of the singular value plateau is crucial. However, 
increasing N at fixed £ max involves only a polynomial 
increase of the computational effort, so we can always 
assume that N is sufficiently large for the plateau to be 
fully developed. 
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FIG. 6: (color online) Comparison for the three-peak spec- 
trum. 

Considering a broader selection of points from G'(r, t) 
data, we did not find an improvement of quality of 
the spectral function compared to pure Keldysh Green's 
functions. 



D. Comparison to direct inversion 

The maximum entropy deconvolution is superior to di- 
rect inversion of the Fourier transform for finite Keldysh 
branches. To demonstrate this, we apply the usual rota- 
tion in Keldysh space and consider the retarded Green's 
function defined in Eq. In contrast to Eq. Q, no 
Fermi factor appears in the transform from the spectral 
function here, 

/oo 
due-W-^Afa), t>t'. (14) 
-oo 

In the limit of a Keldysh contour of infinite length, a 
Laplace transform restores the spectral function 

i r°° 

A(u) = -hn- dfc +ltJt G rct (i). (15) 

7f JO 

For our contour of finite lenght t max a straightforward 
approximation to the deconvolution problem involves a 
truncation of the Laplace integral at time t max . Figure 
[7] shows spectra resulting from this method for the rect- 
angular test spectrum and contour lengths i max = 2 and 

^max 10. 

Comparison to the finite-e maximum entropy results 
from Keldysh propagators in Fig. [4] shows that the sharp 
edge is not resolved as well. Furthermore, unphysical re- 
gions of negative spectral weight appear in the case of 
the truncated Laplace transform, whereas this is avoided 
by construction in the maximum entropy approach. The 
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FIG. 7: (color online) Spectra obtained from retarded Green's 
function on the Keldysh contour using a truncated Laplace 
transform (£ max = 2 and t max = 10) and comparison to the 
exact result. The i max = 10 Maximum Entropy result is in- 
cluded for comparison. See Fig. [4] for a more detailed com- 
parison with Maximum Entropy data. 



amplitude of the ringing oscillations seems to be almost 
the same as in the maximum entropy case, except near 
the band-edge, where the truncation involves an arbi- 
trary smoothening of sharp structures. Different cut- 
off procedures for the time integral would yield different 
smoothenings of the spectrum. 



III. PADE APPROXIMATION 

As a second alternative to the Maximum Entropy ap- 
proach, we consider the Pade method, which in equilib- 
rium situations is often superior to the Maximum En- 
tropy method for the analytic continuation of data with- 
out, or with very little stochastic noise. In particular for 
the analytical continuation of low-temperature Matsub- 
ara data, it is known to be rather precise at low frequen- 
cies [H]. 

The Pade approximant is constructed as follows: We 
assume that the values of the Green's function are known 
on a set of N points z n in the complex frequency plane. 
Conventionally, the z n are given by the Matsubara fre- 
quencies iuj n (uj n = (2n + 1)tv//3). The Green function is 
interpolated with a rational function Cj\r(z) in the form 
of a continued fraction, 

a 1 a 2 (z-z 1 ) a N (z-z N -i) 

<^N(Z) — 7-. r-. ■■■ : , (.J-DJ 



1- 



1- 



1 



where the coefficients of the continued fraction are 
computed using a simple recursion formula [25] . 

In the present paper we are in particular interested 
in the non-equilibrium situation, i.e. the determination 
of a time-dependent spectral function A(t, tu) from the 



FIG. 8: (color online) Approximation of the rectangular spec- 
trum with the Pade approach, either exact or for £ max = 10. 
See Figs. [4] and [7] for comparison with Maximum Entropy and 
the truncated Laplace transform. 



real-time retarded Green's function on some finite time 
interval [t, t + At]. For a set of points z n in the complex 
plane, we first compute an approximation G(t, z n ) of the 
Fourier transformed G Ict , 



G(t, z n 



t+At 



G rel! (t / ,i)e i *» (t '- t) dt' . 



(17) 



This approximation is good if the imaginary part of z n 
is sufficiently large (Imz n ^> l/At). One can then con- 
struct a Pade approximant to the function G(t, z) which 
tends to suppress the artificial oscillatory behaviour of 
spectral functions obtained from direct Fourier trans- 
forms (corresponding to z — > lj + i0 + ). 

For the approximate Pade method we consider the fol- 
lowing three variants: the points z n are 

• "imag" : fermionic Matsubara frequencies z n = 
ia; n =i(2n + l)n/f3, with some arbitrarily chosen /3, 

• "real" : real frequencies with some imaginary offset 
vy, i.e. z n = nAco + vy, 

• "grid": arranged on a N lca \ x JV imag grid, where 



,.(r) 



liO 



(') 



The lattice is chosen 



(r) 

to be equidistant, with uin IB£ 



(i) 



w„ _i = Auj and 



7- 



Figure |8] shows spectra obtained from these different 
Pade variants, when applied to the rectangular test spec- 
trum. For the approximate Pade versions we always use 
the Keldysh contour length t max = 10, as well as 100 
data points as input for any of the Pade procedures. It 
is interesting to note that, as compared to the truncated 
Laplace transform and the Maximum Entropy approach, 
the amplitude of the oscillations is significantly smaller 



for the conventional Pade. However, the resolution of the 
jump is clearly limited. As /? is decreased for the Mat- 
subara frequencies, a slight increase of the slope can be 
observed. 

Let us now discuss the approximate Pade method for 
a time interval length t max — 10. By construction it is 
closely related to the truncated Laplace transform, whose 
results are shown in Fig. [7J Nevertheless, features of the 
original Pade approach are also inherited, depending on 
how the z n are chosen. In fact, the "imag" solution is 
much closer to the conventional Pade solution than to 
the respective truncated Laplace transform. The unde- 
sirable oscillations have been smoothened at the price of 
a reduced resolution at higher frequencies, while some 
of the oscillations from the truncated Laplace transform 
survive. For 7 = 1, the "real" and "grid" solutions are 
almost identical. These results are closer to the exact so- 
lution than the truncated Laplace transform, even near 
the jump, and also exhibit smaller amplitude oscillations. 
Unfortunately, as shown in the lower panels, the rather 
unsystematic behaviour as a function of 7 makes it diffi- 
cult to determine the optimal value of 7 a priori. We also 
note that in the "real" variant, the truncated Laplace so- 
lution is recovered for 7 = 0.5, since in this case all the 
z n are close to the real axis. 



IV. APPLICATION TO DYNAMICAL MEAN 
FIELD RESULTS 

A. Equilibrium spectra 

As a realistic application, we analyze data obtained for 
the Kondo lattice model (bandwidth 4, (3 = 50, J = 1.5, 
3.0) within dynamical mean field theory in Ref. |26j . us- 
ing the non-crossing approximation (NCA) as impurity 
solver. We do not want to discuss here the physics of the 
Kondo lattice model, and to what extent certain high en- 
ergy features in the spectral function may be an artifact 
of the NCA. We merely use the results of Ref. [33] as 
a nontrivial example involving low-energy quasi-particle 
peaks and high-energy satellites. 

In order to apply the Maximum Entropy Method 
(MEM), we again assume a uniform error e and set off- 
diagonal entries of the covariance to zero. This seems 
justified as long as the systematic numerical error of the 
simulation data is significantly smaller than e. In prac- 
tice, the value of e can be rather easily determined by 
trial and error: if e is chosen too small, the MEM ceases 
to converge, whereas if it is chosen too large, the data lack 
sufficient information. In the following calculations, e is 
chosen slightly larger than the value at which the MEM 
breaks down. This moves the systematic numerical errors 
just into the e-range of the Gaussian statistical errors 
which arc formally assumed by the Maximum Entropy 
procedure. As input data we consider the Matsubara 
Green's functions G m (t) and the retarded (equilibrium) 
Green's functions G ret (t — t'), as well as data sets con- 
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FIG. 9: (color online) Comparison of Pade and conventional 
MEM spectra for a Kondo lattice model (bandwidth 4, inverse 
temperature f3 — 50), in the Fermi-liquid regime (J = 1.5, 
left) and Kondo insulating regime (J = 3.0, right). 

taining information from both Matsubara and retarded 
Green's functions. 

The NCA data considered here are noise-free low tem- 
perature data, and thus ideally suited for the Pade 
method. In Fig. [9] we compare MEM, Pade, and ex- 
act spectra for a half-filled Kondo- insulator, and a doped 
heavy-Fermi liquid solution. The "exact" spectra have 
been obtained from a truncated Laplace transformation 
with large i max ~ 45. We see that both Pade and Mat- 
subara MEM reproduce the low-energy features and the 
gap-edges very well, but the higher energy structures can- 
not be accurately resolved. 

In Fig. [lO] we show in addition to the exact and Mat- 
subara MEM spectra the results from a MEM analytical 
continuation of real-time data with i max = 10 and of 
a MEM continuation involving both the real-time and 
the imaginary-time data. The real-time MEM spec- 
tra produce the correct high-frequency behavior, and at 
least qualitatively correct structures in the lower band of 
the doped system, but they fail to reproduce the quasi- 
particle peaks. By adding also the imaginary-time data, 
we can accurately resolve the low-energy behavior (see in- 
sets), in addition to the high-energy structures. The main 
challenge in the equilibrium case thus remains the cal- 
culation of spectral features in the intermediate-energy 
range, where the convergence to the exact result with 
increasing £ max is relatively slow. 



B. Time-dependent spectra 

To test the ability of the MEM approach to extract 
time-dependent spectra, we again consider results from 
Ref. [551 These spectra correspond to a perturbed, doped 
Kondo-lattice model which dissipates energy to a heat- 
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FIG. 10: (color online) Comparison of MEM solutions for 
the insulating Kondo lattice model (top, J = 1.5 and 3.0) 
and the doped heavy Fermi liquid regime (bottom, J = 1.5). 
A real-time contour length t max = 10 is sufficient to extract 
the relevant features. 



bath, and thereby evolves from the high-temperature 
local moment regime into the low-temperature heavy- 
Fermi liquid regime. Again, we are not interested in the 
physics here, which has been discussed in Ref. [26], but 
just in the quality of the spectral functions that can be 
extracted from the retarded Green functions. In partic- 
ular, we want to investigate the constraints on t max for 
the extraction of a time-dependent spectral function at 
time t {t < i m ax), since most computational techniques 
are limited to short time contours. 

The upper panels of Fig. [IT] show the spectral func- 
tion A(io, t) at a relatively short time t = 6 after the 
perturbation. As can be seen, A(ui,t) is already positive 
over the whole relevant frequency range. The different 
curves in the figure correspond to different time windows 
[t, t + At] used in the MEM analytical continuation. The 
result for the largest At (dash-dotted curve) can be con- 
sidered the exact result. We learn from these results that 
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FIG. 11: (color online) Convergence of A(tJ, t — 6) for differ- 
ent time-intervals [t, t + At]. 



the high-frequency part converges very quickly with At, 
while we need at least an interval of length At — 24 to 
get a reasonably accurate result at low energies. (Since 
this calculation is a real nonequilibrium application, we 
cannot resort to Matsubara data to fix the low-energy 
part of the spectral function.) 

A relevant question is how the MEM approach per- 
forms compared to direct Fourier transformation on an 
identical time interval. While both the MEM method and 
the direct Fourier transformation have difficulties repro- 
ducing the correct low-energy peak, the MEM spectra 
are slightly closer to the correct result. 



C. Convergence of high-energy features 

While the convergence of the MEM result for pure real- 
time Green's function data as a function of At is slow 
for the quasi-particle peak, it is interesting to see how 
the higher-energy features of the non-equilibrium spec- 
tral function converge to the exact result. 



Figure 12 shows a comparison of the MEM result to the 
direct Fourier transform and the modified Pade method 



(discussed in Sec. Ill ) for short time intervals with lengths 
At e [3.9,5.7]. For At = 5.7, MEM and direct Fourier 
transformation produce the correct solution, except in 
the low-energy region —0.5 < u> < 0.5, the modified Pade 
method for the imaginary time interval yields spurious 
results. For shorter interval lengths, the MEM appears 
to be slightly superior to both, direct Fourier transform 
and modified Pade. All qualitative features of the spec- 
trum, including the double-peak structure in the lower 
band start to become visible in the MEM solution at 
At = 4.5, whereas the emergence of similar structures in 
the direct FT spectrum at At = 5.1 may still be inter- 
preted as part of the overall oscillatory behaviour of that 
solution. One may thus argue that the MEM reduces 
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FIG. 12: Convergence of non-equilibrium A(uj,t = 6) estimates for different short real-time intervals. In the case of the "imag" 
modified Pade, we set ft = 20 and used 100 data points. For the "real" modified Pade results, we used 100 equidistant data 
points for Rez n on the interval [—10, 10]. 



the necessary interval length from At = 5.7 to At = 4.5 
and thus allows to compute time-dependent spectra up 
so slightly longer times. Depending on the application, 
this could safe considerable effort, since numerical algo- 
rithms typically scale badly (at least as the third power) 
with the length of the contour. The effect of the modified 
Pade for z n on parallel lines to the real axis can to a large 
extent be interpreted as a smoothening procedure for the 
direct FT. In particular, the double-peak structure in the 
direct FT solution at At = 5.1 is strongly suppressed, 
but still slightly visible. However, the peaks are less pro- 
nounced than those of the MEM solution at At = 4.5. 
The behaviour of the "grid" variant is very similar to the 
"real" variant, except that the "grid" choice cannot re- 
solve the double-peak structure of the lower band (not 
shown) . When the number of real frequency points is in- 
creased, the "grid" result converges to the "real" result. 
Due to numerical instabilities of the Pade procedure, one 
is however typically limited to less than 400 grid points. 



V. CONCLUSION 

We analyzed the usefulness of Kcldysh real-time 
Green's function data for the computation of equilib- 
rium spectral functions of interacting quantum many- 
body systems. The ill-posed nature of the inversion cor- 
responding to the determination of the spectral function 
is reflected in the distribution of singular values of the 
respective continuation kernel. The conventional Wick 
rotation of Matsubara data yields an exponentially de- 
caying singular value distribution. We found that in- 
cluding a finite real-time branch to the imaginary-time 
contour adds a plateau to this distribution, which broad- 
ens as the length of the real-time contour is increased. 
This plateau significantly alleviates the inversion even 



for small contour lengths. A further analysis showed that 
the main gain is more precise high-frequency information 
of the spectral function. For short contours, the real- 
time data however provide rather crude results at low 
frequencies. In order to obtain more accurate estimates 
of low-energy features such as quasi-particle peaks, it is 
therefore necessary to include data from the Matsubara 
branch. With the small frequencies well covered by the 
Matsubara branch and the high frequencies well covered 
by a short Keldysh branch, the biggest challenge remains 
the resolution at intermediate energies. 

Our analysis of NCA spectra for the Kondo-lattice 
model suggests that in order to resolve intermediate en- 
ergies reasonably well, the length of the real-time con- 
tour has to be at least t max — 10 (in units where the 
bandwidth is 4). In real-time Monte Carlo simulations, 
it is difficult to reach these times, at least in parame- 
ter regimes where low-order perturbation theory is not 
reliable. It thus appears that the exponential increase 
in the computational cost with increasing £ max in real- 
time Monte Carlo schemes outweighs the benefits from 
an alleviated MEM analytical continuation. Neverthe- 
less, the MEM analytical continuation of Keldysh Green 
functions can be a useful and superior alternative to Pade 
approximants for semi-analytic methods such as higher 
order strong-coupling perturbation theory. These meth- 
ods are very useful in certain parameter regimes, such as 
elevated temperature, or intermediate to strong interac- 
tions. While their implementation on the real-frequency 
axis is challenging, calculations on the Keldysh contour 
scale polynomially with i max and the contour-lengths 
needed for reliable MEM analytical continuation can be 
reached. 

In the calculation of non-equilibrium spectra, it is no 
longer possible to include the Matsubara branch into the 
Maximum Entropy procedure. As a consequence, the 
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resolution of low-energy features of the spectrum wors- 
ens dramatically if only short-time data are available. 
Nevertheless, a comparison to alternative, more straight- 
forward techniques, i.e. direct Laplace transform and the 
modified Pade approach, showed that the MEM yields 
slightly more reliable solutions for these spectra. The 
relevant structures could be established for somewhat 
shorter time intervals than for the extended Pade ap- 
proaches. 

On a conceptual level, the MEM is superior to the 
direct Laplace transform and the generalized Pade ap- 
proach. In practical applications, however, because the 
advantage of a MEM continuation appears to be rather 
subtle, it will very much depend on the details of the 
utilized many-body approach whether an application of 



MEM is worth its effort. In any case, the comparison of 
different approaches can help in deciding which features 
of a non-equilibrium spectrum are trustworthy. 
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